Zadanie 1

Porównać w języku Julia reprezentację bitową liczby 1/3 dla Float16, Float32, Float64 oraz liczby, która jest inicjalizowana jako Float16, a potem rzutowana na Float64.

Rozwiązanie

#using Pkg
#Pkg.add("DataFrames")
using DataFrames

fl = [Float16(1/3), Float32(1/3), Float64(1/3), Float64(Float16(1/3))]
fl_bit = map(x -> bitstring(x), fl)

floats = DataFrame()
floats.Float_type=["Float16", "Float32", "Float64", "Float16->Float64"]
floats.Sign_bit = map(x -> x[1], fl_bit)
floats.Exponent_bits = [fl_bit[1][2:6], fl_bit[2][2:9], fl_bit[3][2:12], fl_bit[4][2:12]]
floats.Mantissa_bits = [fl_bit[1][7:16], fl_bit[2][10:32], fl_bit[3][13:64], fl_bit[4][13:64]]

floats

4×4 DataFrame

Row	Float_type	Sign_bit	Exponent_bits	Mantissa_bits
	String	Char	String	String
1	Float16	0	01111	1111010101
2	Float32	0	01111111	10101010101010101010101
3	Float64	0	01111111101	0101010101010101010101010101010101010101010101010101
4	Float16->Float64	0	01111111101	0101010101000000000000000000000000000000000000000000

Okazuje się, że po rzutowaniu Float16 na Float64 pozostała część mantysy wypełniana jest zerami. Naturalnie taka liczba mniej dokładnie oddaje wartość 1/3 niż liczba pierwotnie zainicjalizowana jako Float64.

Zadanie 2

Zbadać, jak zmienia się odległość między kolejnymi liczbami zminnoprzecinkowymi reprezentowanymi w komputerze za pomocą języka Julia. Narysować wykres używając Plots zależności odległości od wartości liczby dla zakresu od 1.0 do 1000000.0.

Rozwiązanie

#using Pkg
#Pkg.add("Plots")
using Plots
r = 1:1000000
epsilons = [eps(Float64(n)) for n in r]
plot(r, epsilons, title="", label="", xlabel="value of the number", ylabel="distance between neighbouring numbers")

Zadanie 3

Jedną z bibliotek numerycznych, jaką dodatkowo będziemy używać na zajęciach jest GSL (język C). Opis jak używać . Korzystając ze wsparcia dla wyświetlania reprezentacji liczb zmiennoprzecinkowych zobaczyć jak zmienia się cecha i mantysa dla coraz mniejszych liczb. Zaobserwować, kiedy matysa przestaje być znormalizowana i dlaczego?

Kod załączyć jako komórka Markdown sformatowana jako C (link). Wynik także jako Markdown (kod albo fragment zrzutu ekranu).

Rozwiązanie

Kod w C z użyciem biblioteki gls:

#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>
#include <math.h>

int main()
{
    // floats are 32bit
    float f = 1.12345678 * exp2(-124);

    for (size_t i = 0; i < 27; i++)
    {
    	f/=2;
    	gsl_ieee_printf_float(&f);
    	printf("\n");
    }
    
    return 0;
}

Poniżej wydruk konsoli:

1.00011111100110101101111*2^-125
 1.00011111100110101101111*2^-126
 0.10001111110011010111000*2^-126
 0.01000111111001101011100*2^-126
 0.00100011111100110101110*2^-126
 0.00010001111110011010111*2^-126
 0.00001000111111001101100*2^-126
 0.00000100011111100110110*2^-126
 0.00000010001111110011011*2^-126
 0.00000001000111111001110*2^-126
 0.00000000100011111100111*2^-126
 0.00000000010001111110100*2^-126
 0.00000000001000111111010*2^-126
 0.00000000000100011111101*2^-126
 0.00000000000010001111110*2^-126
 0.00000000000001000111111*2^-126
 0.00000000000000100100000*2^-126
 0.00000000000000010010000*2^-126
 0.00000000000000001001000*2^-126
 0.00000000000000000100100*2^-126
 0.00000000000000000010010*2^-126
 0.00000000000000000001001*2^-126
 0.00000000000000000000100*2^-126
 0.00000000000000000000010*2^-126
 0.00000000000000000000001*2^-126
 0
 0

0

Mantysa przestaje być znormalizowana, gdy liczba, którą system musi przechować jest mniejsza niż $2^{-126} \cdot 1.0$. Wynika to z tego, że minimalna wartość mantysy to $1$ (dzięki domyślnej jedynce), a najmniejsza wartość cechy bez wszystkich zerowych bitów (aby była znormalizowana) to 1 (00000001), a po uwzględnieniu $biasu=-127$ ostatecznie mamy $2^{1-127}=2^{-126}$. Poniżej tej wartości granicznej liczby muszą być reprezentowane w postaci zdenormalizowanej (cecha ma postać (00000000)) a mantysa traci domyślną jedynkę, przez co należy teraz do przedziału (0, 1). Najmniejszą możliwą do przechowania liczbą w takim układzie jest $(0,00000000, 00000000000000000000001)$

Zadanie 4

Na przykładzie wybranego algorytmu niestabilnego numerycznie:

1. Pokazać, że działa źle.
2. Pokazać które konkretnie działania powodują zwiększenie błędu (np. dzielenie przez małą liczbę, cancellation).
3. Zademonstować wersję stabilną.

Wszystkie punkty przedstawić w postaci notatnika Julii.

Rozwiązanie

Jako nasze zadanie przyjmijmy obliczenie wartości całek $$y_n = \int_{0}^{1} \frac{x^n}{x+5} dx$$ dla n=0, 1, 2, ... 15.

Zauważmy, że $$ y_n + 5y_{n-1} = \int_0^1 \frac{x^n + 5x^{n-1}}{x+5} dx = \int_0^1 \frac{x^{n-1}(x+5)}{x+5} dx = \left [ \frac{x^n}{n} \right ] _0 ^1 = \frac{1}{n} $$ co daje nam wzór rekurencyjny $$y_n + 5y_{n-1} = \frac{1}{n}$$ Po przekształceniu wzoru otrzymujemy zależność $$y_n = \frac{1}{n} - 5y_{n-1}$$ Potrzebujemy jeszcze wartości pierwszego wyrazu $$y_0 = \int_0^1 \frac{1}{x+5} = \left [ \ln (x+5) \right ] _0^1 = \ln 6 - \ln 5 \approx 0.18232156$$

Przechodzimy do obliczania kolejnych $y_n$

1)

numOfSteps = 15
y = Vector{Float32}()
append!(y, 0.18232156)

for n in 1:numOfSteps
    append!(y, 1/n - 5*last(y))
end
results = DataFrame(y_index = [i for i in 0:numOfSteps], value = y)

16×2 DataFrame

Row	y_index	value
	Int64	Float32
1	0	0.182322
2	1	0.0883922
3	2	0.058039
4	3	0.0431383
5	4	0.0343086
6	5	0.0284571
7	6	0.0243812
8	7	0.020951
9	8	0.020245
10	9	0.00988603
11	10	0.0505698
12	11	-0.16194
13	12	0.893034
14	13	-4.38825
15	14	22.0127
16	15	-109.997

Dla $y_{11}$ otrzymujemy oczywistą sprzeczność, gdyż całka oznaczona w granicach od 0 do 1 z funkcji nieujemnej $\frac{x^n}{x+5}$ nie może mieć wartości mniejszej od zera

2)

Powyższy algorytm z każdą iteracją zwiększa błąd, wynikający z niedosknałości reprezentacji liczb zmiennoprzecinkowych w komputerze. Winę za to ponosi wykorzystana przez nas zależność rekurencyjna $$y_n = \frac{1}{n} - 5y_{n-1}$$ Wynika z niej, że z każdą iteracją mnożymy już posiadany błąd razy 5, więc, w szczególności, pierwszy popełniony błąd $\epsilon_0$ (podczas zaokrąglenia $y_0$) będzie już znacząco wpływać na wartość np. wyrażenia $y_8$, bo będzie to $5^8 \cdot \epsilon_0$!

3)

Aby algorytm był jednak stabilny, zastosujemy inne podejście. Przekształcamy znaną już zależność rekurencyjną do postaci $$y_{n-1} = \frac{1}{5n} - \frac{1}{5} y_n$$

Teraz błąd będzie z każdym krokiem dzielony przez -5. Ponieważ $y_n$ maleje gdy $n$ rośnie, możemy przypuszczać, że dla dużych $n$, $y_n$ maleje wolno, więc przyjmijmy, że $y_{16} \approx y_{17}$ i korzystając z wzoru $$y_{16} = \frac{1}{5 \cdot 17} - \frac{1}{5}y_{17}$$ otrzymujemy $$y_{16} \approx \frac{1}{5 \cdot 17} - \frac{1}{5} y_{16} \Rightarrow y_{16} \approx \frac{1}{5 \cdot 17} \cdot \frac{5}{6} \approx 0.009803921$$

Teraz stosujemy analogiczne podejście, jak w poprzednim algorytmie, ale zaczynając oczywiście od $y_{16}$ a na $y_0$ kończąc

numOfSteps = 16
y = Vector{Float32}()
append!(y, 0.009803921)

for n in numOfSteps:-1:1
    append!(y, 1/(5*n) - last(y)/5)
end
results = DataFrame(y_index = [i for i in numOfSteps:-1:0], value = y)

17×2 DataFrame

Row	y_index	value
	Int64	Float32
1	16	0.00980392
2	15	0.0105392
3	14	0.0112255
4	13	0.0120406
5	12	0.0129765
6	11	0.0140714
7	10	0.0153675
8	9	0.0169265
9	8	0.0188369
10	7	0.0212326
11	6	0.0243249
12	5	0.0284684
13	4	0.0343063
14	3	0.0431387
15	2	0.0580389
16	1	0.0883922
17	0	0.182322

Tym razem wartość $y_0 \approx 0.182322$, co jest poprawnym wynikiem